% current_tfp.m
%
% extracts current TFP shock from a Cholesky decomposition
% i.e. implicit assumption is that only two shocks can affect TFP and that
% news shock has no contemporaneous effect on TFP
%
% Andre Kurmann, Federal Reserve Board; last modified December 2012
%-------------------------------------------------------------------------


function [output,v]=RecursivenessIdentification(b,res,vmat,nvars,nlags,nimp,pos_shock)
% PURPOSE: extracts shock from a Cholesky decomposition; (e.g. current TFP, pos_shock = 1)
% -----------------------------------------------------
% USAGE: [output,v]=RecursivenessIdentification(b,res,vmat,nvars,nlags,nimp,pos_shock)
% where: 
% b    = (3 x 1) vector with indexes for VAR specification: nm,nd,nlend 
% res  = (ntall x (nm+nd) data matrix = [Y D]  
% vmat
% nvars
% nlags
% nimp
% pos_shock
% NOTE: no dimension checks imposed , 31/01/04
% -----------------------------------------------------
% RETURNS: 
% output = (nt x nm) matrix with obss on endog variables
% v      = (nt x nkk) matrix with regressors
% 
% -----------------------------------------------------

%companion matrix of demeaned VAR (=> don't include constant term)
M=zeros(nvars*nlags,nvars*nlags);			
M(1:nvars,:)=b(1:nvars*nlags,:)';
M(nvars+1:nvars*nlags,1:nvars*nlags-nvars)=eye(nvars*nlags-nvars);

% Extract impulse vectors 
%------------------------

%unique lower triangular matrix = Cholesky decomp of vmat
  Gtilde =                              chol(vmat)'; % Impulse Responses to a One Standard Deviation
% Gtilde = (inv(diag(diag(chol(vmat))))*chol(vmat))';% Impulse Responses to a One Percent Innovation

gamma = zeros(nvars,1); gamma(pos_shock) = 1;  %'selection vector' of same format than in sims extraction
                                               %#ATTENTION: Ramey News must be ordered "pos_shock"
alpha = Gtilde*gamma;                          %extracts first column of cholesky

% Back out time series of structural shock 
v = gamma'*inv(Gtilde)*res';    %1xT vector of structural (standardized) shocks
v = v';                         %-->      gamma'*Gtilde      gives the volatility of the structural shock
                                %--> sqrt(gamma'*vmat*gamma)
    
    % Compute impulse responses
    %----------------------------------------------------------------------
    U1=[alpha; zeros(nvars*nlags-nvars,1)];
%     n = size(U1,1); 
%     nlong = 20;   %ATTENTION: specify number of quarters for long rate: for 5-year = 20

%     if use_yields == 0 | use_yields == 3;
        for k=1:nimp
            Zk1(k,:)=(M^(k-1)*U1)';
        end
        
%     elseif use_yields == 1
%         hn = zeros(1,n); h1 = zeros(1,n);
%         if slope == 0;
%             hn(nvars-1)=1;       %selection vector for long rate; long rate is second-to-last of nvars elements if slope=0
%             h1(nvars-1)=1; h1(nvars)=-1;     %short rate = long rate - spread; spread is last of nvars elements if slope=0
%         elseif slope == 1;
%             hn(nvars)=1;       %selection vector for long rate; long rate is last of nvars elements if slope=1
%             h1(nvars)=1; h1(nvars-1)=-1;     %short rate = long rate - spread; spread is second to last of nvars elements if slope=1
%         end
%         hinf = zeros(1,n); hinf(nvars-3)=1;   %inflation; is usually fourth to last in this specification; ATTENTION: ONLY WORKS IF INFLATION IS IN APPROPRIATE POSITIION IN SIMS-VAR!!
%     
%         for k=1:nimp;
%             Zk1(k,:)=(M^(k-1)*U1)';
%             LongEH(k,1) = 1/nlong*(h1*inv(eye(n)-M)*(eye(n)-M^nlong)*Zk1(k,:)')'; %long-rate reaction to 1st shock according to Expectations Hypothesis
%             TP(k,1) = hn*Zk1(k,:)' - LongEH(k); %term-premium reaction to 1st shock
%             SpreadEH(k,1) = LongEH(k,1) - h1*Zk1(k,:)';
%         end
%         
%     elseif use_yields == 2
%             hlong = zeros(1,n); hlong(nvars) = 1;   %long bond yield is last element in sims-var if slope=0 
%             hspread = zeros(1,n); hspread(nvars) = 1; %long bond yield is last element in sims-var if slope=1
%             hffr = zeros(1,n);  hffr(nvars-1) = 1;    %selection vector for short rate = ffr; ATTENTION: ONLY WORKS IF FFR IS IN SECOND-TO-LAST POSITIION IN SIMS-VAR!!
%             hinf = zeros(1,n); hinf(nvars-2)=1;   %selection vector for inflation; ATTENTION: ONLY WORKS IF INFLATION IS IN THIRD-TO-LAST POSITIION IN SIMS-VAR!!
%             for k=1:nimp;
%                 Zk1(k,:)=(M^(k-1)*U1)';
%                 if slope == 0;
%                     Spread(k,1)=(hlong-hffr)*Zk1(k,:)';     %Spread computed from long - ffr
%                     Long(k,1)=hlong*Zk1(k,:)';
%                 elseif slope == 1;
%                     Long(k,1) = (hspread+hffr)*Zk1(k,:)';   %Long rate computed from spread + ffr
%                 end
%                 LongEH(k,1) = 1/nlong*(hffr*inv(eye(n)-M)*(eye(n)-M^nlong)*Zk1(k,:)')'; %long-rate reaction to 1st shock according to Expectations Hypothesis
%                 TP(k,1) = Long(k,1) - LongEH(k);    %term-premium reaction to 1st shock
%                 SpreadEH(k,1) = LongEH(k,1) - hffr*Zk1(k,:)';
%             end
%     end
    
    impulse1(:,:)=Zk1(:,1:nvars);
%     if use_yields == 1 | use_yields == 2
%        impulse_add1(:,1) = LongEH;  
%        impulse_add1(:,2) = TP;
%        impulse_add1(:,3) = SpreadEH;
%        if use_yields == 2
%             if slope == 0
%                    impulse_add1(:,4) = Spread;
%             elseif slope == 1;
%                    impulse_add1(:,4) = Long;
%             end
%        end
%     end 
    
    %test whether TFP impulse response is positive on impact, otherwise revert
    %#ATTENTION: NEED TO CORRECTLY SPECIFY COLUMN OF RELEVANT VARIABLE !!!
    %TFP is target for FEV maximization and TFP is 1st in VAR
    if Zk1(1,pos_shock) < 0   
       alpha = - alpha; 
       v = -v;
       impulse1 = -impulse1;
%        if use_yields == 1 | use_yields == 2;
%             impulse_add1 = -impulse_add1;
%        end
    end

    
    % Compute fraction of VD explained at different horizons
    %----------------------------------------------------------------------
    B = zeros(nvars,nvars,nimp+1);
    for l=0:nimp
        C=M^l;
        B(:,:,l+1) = C(1:nvars,1:nvars);
    end
    sigmak=B(:,:,1)*vmat*B(:,:,1)';
    hh1=B(:,:,1)*alpha*(B(:,:,1)*alpha)';
    vardec1(1,:)=(diag(hh1./sigmak))';
%     if use_yields == 2;
%         if slope == 0;
%             hlong = zeros(1,nvars); hlong(nvars) = 1;   %long bond yield is last element in Sims-VAR
%             hffr = zeros(1,nvars);  hffr(nvars-1) = 1;    %selection vector for short rate = ffr; ATTENTION: ONLY WORKS IF FFR IS IN SECOND-TO-LAST POSITIION IN SIMS-VAR!!
%             sigmak_spread = (hlong-hffr)*sigmak*(hlong-hffr)';    %FEV of spread computed from FEV(long-ffr)
%             hh1_spread = (hlong-hffr)*hh1*(hlong-hffr)';
%             vardec1_add(1,1) = hh1_spread/sigmak_spread; 
%         elseif slope == 1;
%             hspread = zeros(1,nvars); hspread(nvars) = 1;   %long bond yield is last element in Sims-VAR 
%             hffr = zeros(1,nvars);  hffr(nvars-1) = 1;    %selection vector for short rate = ffr; ATTENTION: ONLY WORKS IF FFR IS IN SECOND-TO-LAST POSITIION IN SIMS-VAR!!          
%             sigmak_long = (hspread+hffr)*sigmak*(hspread+hffr)';    %FEV of long rate computed from FEV(spread+ffr)
%             hh1_long = (hspread+hffr)*hh1*(hspread+hffr)';
%             vardec1_add(1,1) = hh1_long/sigmak_long;
%         end
%     end

    for k=1:nimp-1
        %add square of k-step ahead forecast error to build k-ahead variance-covariance (eq. 6)
        sigmak=sigmak+B(:,:,k+1)*vmat*B(:,:,k+1)';
        %fraction explained by the news shock (eq 8)    % XX change in descripion (AK)
        hh1=hh1+B(:,:,k+1)*alpha*(B(:,:,k+1)*alpha)';
        vardec1(k+1,:)=(diag(hh1./sigmak))';
%         if use_yields == 2;
%             if slope == 0;
%                 sigmak_spread = (hlong-hffr)*sigmak*(hlong-hffr)';    %FEV of spread computed from FEV(long-ffr)
%                 hh1_spread = (hlong-hffr)*hh1*(hlong-hffr)';
%                 vardec1_add(k+1,1) = hh1_spread/sigmak_spread;    
%             elseif slope == 1;
%                 sigmak_long = (hspread+hffr)*sigmak*(hspread+hffr)';    %FEV of long rate computed from FEV(spread+ffr)
%                 hh1_long = (hspread+hffr)*hh1*(hspread+hffr)';
%                 vardec1_add(k+1,1) = hh1_long/sigmak_long;
%             end
%         end
    end      
    

    output=[vardec1 impulse1];
